STATISTICAL EXTRAPOLATION 



Dennis Randolph Oldson 



United States 
Naval Postgraduate School 




STATISTICAL EXTRAPOLATION 




by 


Dennis 


Randolph Oldson 


Thesis Advisor: 


Sydney R. Parker 



September 1971 



Apptoue.d puhtlc ■tc.le.ase.; dUi.'Ubution mUjnite.d. 



/Aes/'s 



Statistical Extrapolation 



by 



Dennis Randolph Oldson 
Lieutenant, United States Navy 
B.S.E.E., University of Colorado, 1965 



Submitted in partial fulfillment of the 
requirements for the degree of 



ELECTRICAL ENGINEER 



from the 

NAVAL POSTGRADUATE SCHOOL 
September 1971 



ABSTRACT 



The mathematics of extrapolating known statistics of 
components to the probability density function of a system's 
performance measure is considered. Quadrature sum inte- 
gration schemes for evaluating the resulting required inte- 
gration are examined, and alternate integral approximation 
schemes are developed utilizing Monte Carlo methods. A 
simple electrical circuit example illustrates the use of 
these techniques. 
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I. 



INTRODUCTION 



The desirability of being able to statistically predict 
a system's performance measure is obvious when one considers 
the mass production of systems. Since performance measures 
of certain values may not be desirable, it is generally un- 
profitable to produce such systems when the probability that 
these undesirable values can exist is beyond the limits 
determined by the requirements of the situation. 

Presently performance data for systems are obtained in 
several ways: (1) A worst case analysis [1] to determine if 

undesirable performance will occur within the parameter toler- 
ance; (2) A moment method [1] in which the performance 
measure is sampled so that its mean, variance, and higher 
order moments are determined for approximation by known dis- 
tribution functions; (3) A functional formulas method [2] in 
which performance measure's distributions are obtained from 
breaking performance measures into a series of products and 
sums for individual integration; (4) Monte Carlo methods 
[3, 10] in which performance measure distributions are gener- 
ated by creation of random processes by which the parameters 
of the random process lead to an approximation to the desired 
distribution . 

The purpose of this thesis is to investigate the mathe- 
matics required in obtaining a performance measure and at- 
tempt to arrive at a reasonably fast and accurate method of 
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extrapolating the statistics of a system's component variables 
to create the distribution of the desired performance measure. 

Chapter II examines the mathematical form and statistical 
nature of the performance measure for implicit and explicit 
functions and presents two digital computer integration 
schemes using Gauss quadrature sums for evaluating the re- 
sulting required integration. 

Chapter III investigates crude Monte Carlo methods for 
determining performance measures, and the confidence limits 
involved in the method. Then Monte Carlo schemes to estimate 
single integrals are examined. The multiple integral exten- 
sions of the single integral estimates are then derived. 

Chapter IV discusses the relative merits of the methods 
outlined in Chapters II and III, and gives a comparison of 
the results of these methods applied to a simple example con- 
sisting of an electrical circuit in which the component 
values are the random variables. 
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II. STATISTICAL ANALYSIS 



The probability density function of a system's perform- 
ance measure can be determined mathematically from knowledge 
of the functional relation between the performance measure 
and the system component-variables, and knowledge of the 
joint probability density function of the component- 
variables. Most often this requires integrating a complex 
function over several variables. Although the integration 
can be performed analytically in some cases, digital computer 
integration schemes must normally be utilized. 

This chapter examines the mathematics of determining the 
probability density function of a performance measure, first 
for performance measures which are explicit functions of the 
component variables, then for implicit performance measures. 
Two integration schemes are then presented for performing the 
integration on a digital computer. 

A. EXPLICIT CASE 

Consider a performance measure, Z as a function of the 

T T . 

system's component -vari ab le s , X , where X is the transpose 

of an n-vector of elements X^, X^t .../ X^ . It is assumed 

that the X^'s are statistically independent, so that the 

T 

joint probability density function of X is given by the 
product of the probability density functions of the individ- 
ual X . ' s . 

1 
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T 

The probability density function of Z = g(X ) can be 
written as [7] 

T ( X z ) 

where p , Z — ' is the joint probability density function 

A 

T 

for X and Z , V is the n-dimensional space determined by the 

component-variables, g is the functional relation between Z 

T 

and the component variables, X , and dV = dx^dX 2 . . . dx^ . 

If a transformation of variables is defined such that 

= (X2/X3, . . . ,X^) ; ^ = (X^Y^) (2) 

X^ = f(Z,Y^) (3) 

where f is the functional relation between the component- 
variable X^ and the performance measure Z, equation (1) 
becomes [7, 9] 

P 2 (z) = I p^T [ f (z ,y^) ,y^] 1 J (z) 1 dU (4) 

where U is the n-1 dimensional space determined by the n-1 
component -variables Y , dU = dy^dy 2 . . • dy^_^ , y^ = + 

J(z) is the Jacobian of the transformation (2) and (3) . 






(x"^, z) dV 



( 1 ) 
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J{z) can be written as 



J(z) 



3 z 



9z 



3x. 






n-1 



3x 

3z 

3x 

Tz 



n 

1 



0 




3x 

9y 



n 

n-1 



0 





(5) 



where I is the n-1 by n-1 identity matrix. 

Since the joint probability density function of statis- 
tically independent random variables is the product of the 
individual probability density functions, equation (4) 
becomes 






p [f (z,^^) ,y^] j3x,/3zl PyT(y'^)dU 
Ju ^1 ^ L 



( 6 ) 



In the case where the performance measure is a known or 
explicit function of the component-variables, the inverse 
functional relation f relating to Z can usually be ex- 
pressed analytically as can the partial derivative of Z with 
respect to or any of the other component-variables. Cir- 
cumstances may require reordering the component-variables 
to obtain these relations in a convenient form. 
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As an example, consider the series R-L-C circuit in 



Figure 1. The current I through the circuit is 



I ( jto) 



E(jco) 

R + j (coL - 1 /ojC) 



( 7 ) 



and is a maximum for toL = 1/toC. The angular frequency at 
which the current falls to 3 db of the maximum occurs when 



R = toL - 1/toC 



( 8 ) 



or, 

“Sdb ^ ^ R^/4L^ + 1/LC (9) 

Let the performance measure Z be the upper 3 db frequency 
and = (R,L,C) , so that 

Z = g(x'^) = X^/2X2 + / + 1/X^X^ (10) 

and from ( 8) , 

X^ = f(Z,y'^) = ZX 2 - l/ZX^ (11) 

so that the partial derivative is 

8X^/3Z = X 2 + l/Z^X^ (12) 

If each of the components has a Gaussian distribution 
about its mean such that 

exp [-(x .-X .) ^/2a .^] 

Px. = 

3 (2it) ' 

where x^ is the expected or average value of the variable 

X., and a. is the standard deviation of the variable X.. 

3 3 3 
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T T 

Substituting for X, in (13) for j = 1/ P„ If(z,y ) ,y ] 

1 " 

becomes 



— 2 2 

exp [- ( zx 2 ~l/zx^-Xj^) /2a^ ] 



1 



X, 



777172“ 

(2tt) a- 



( 14 ) 



Since X_ and X are independent, p T(y ) is the product of 

Z, ^ 2 

p^ (x„) and p (x^) , and (6) becomes 
X 2 2 ^3 ^ 



exp [ 



p^iz) = 



— 2 

- (zX2-l/zx2“X^) 



- 2 
2 a-, 



/ - ^ 2 

(x^-x^) 

^ 2 

2a-, 



— 00 —00 



(2tt) a ^0202 t (x2+l/z^X2) ] 



-dx dx 
2 3 



(15) 



With a-priori knowledge of the probability density func- 
tions of the component-variables, (6) or (15) can be inte- 
grated on a digital computer using one of many algorithms 
available. The two quadrature formulas presented below offer 
routines which are fast, and for the number of points at 
which the integrand is evaluated, two of the most accurate 
routines available [5, 8]. 



B. IMPLICIT CASE 

If the performance measure is not known explicitly, but 
is obtained implicitly from a solution of a mathematical 
model of the system, relation (3) cannot be determined. One 
may, however, resort to the cumulative distribution obtained 
from (1) by [7] 



✓ 
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F(z) 



P[Z < z] 




V 



(x ,z) dVdz 



T 



(16) 



Since 




(17) 




either 0 or 1, the cumulative distribution can be approxi- 
mated by 



where the integration over the space V may be performed 
utilizing the quadrature formulas below with the value of 



The cumulative distribution function for the performance 
measure can then be fitted to a polynomial, and if desired, 
the polynomial can be differentiated analytically to obtain 
the probability density function. 

C. QUADRATURE SUMS FOR EVALUATING INTEGRALS 

The form of the performance measure has been reduced to 
an integral or multiple integral which, in general, cannot be 
evaluated in closed form. All algorithms which approximate 
the value of an integral by a linear combination of values of 
the integrand are exact [5] for the integrand being of a 



F(z) 



V f T 

I Az p T (x ) dV 

Z<z •'V - 



(18) 



T T T 

PyT (x ) as calculated for Z (x )±z and p T(x ) = 0 if 

A — — A — 



Z ( x"^ ) > z . 
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certain degree or less; and in most cases the degree is one 
less then the number of points at which the integrand is 
evaluated. Quadrature sum algorithms utilize orthoganal 
polynomials to arrive at the form of the linear combination 
of values of the integrand; and as a result are exact for the 
integrand being a polynomial of degree one less than twice 
the number of points at which the integrand is evaluated. 

[5] 

Consider the integral 

fh 

I = h(x)g(x)dx (19) 

•’ a 

where a and b are real numbers, finite or infinite, g(x) is 
an arbitrary function and h(x) is a particular weighting 
function to be described. Then by selecting a polynomial 
Q^(x) of degree n such that h(x)Q^(x) is orthoganal to x^, 
for all m = 0, 1, ..., n-1, so that 

/•b 

x"Vi(x)Q^(x) dx = 0 (20) 

•’ a 



If (20) is approximated by a linear combination of the 
integrand 



' x'\(x)Q_,(x)dx = J,VkX<--^k> 

3. iC — 



then the right hand side of (21) will be equal 
any set of values of Aj^ if the ' s are chosen 



are the zeros of Q (x) . 

n 



( 21 ) 

to zero for 
such that they 
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awigiiw 
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SHPP 






J 












By forming the n sums 



/■b n 

x"Vi(x)dx = I A, 
•’a k=l ^ 



X, 



m 



( 22 ) 



for m = and for the n values of Xj^ being the n 

zeros of Q^(x), the n values of can be found such that 



D n 

h(x)g(x)dx = I A, g(x, ) 
a k=l 



(23) 



is exact for all polynomials of degree 2n-l. [5] 

It can be shown that if the polynomial Q^(x) is known 
then the values of Aj^ can be found to be [5, 8]. 



= 2/(l-Xj^^) [dQ^ (Xj^) /dx] ^ 



(24) 



1 . Gauss Legendre Quadrature 

The first obvious choice of w(x) is 1, and if the 

j 

interval [a, b] can be normalized to the interval [-1, 1] by 
a change of variables, (23) becomes 
fl 



n 



-1 



g(x)dx = I Aj^g(Xj^) 



(25) 



k=l 



and the polynomial (x) is a Legendre polynomial of degree n 



[5] . 



Q (x) = 

n 



,n , 2 T > n 
d (x -1) 

dx 2 n 1 



(26) 



1‘he n values of Xj^ can be found from the zeros of 
Q^(x) and the values A^ from (24) or the values can be ob- 
tained from Table II for values of n from 2 to 6 and from 
15] for values of n from 2 to 48- The magnitude of the error 
bound can be shown to be 15], 
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^rr ^ ^ max 

— (2n+l) (2n) 1 x 



df 



2n 



dx 



2n' 



( 27 ) 



2 . Gauss Hermite Quadrature 

As is often the case when dealing with probability 
density functions, a variable is Gaussian in distribution, 
and as a result an integral of the form 

r" 2 ^ 

h(x)exp(-x )dx = Y A,h(x, ) (2 8) 

J k=l ^ ^ 

— OO 

must be evaluated. This is the form of (19) where w(x) = 
exp ( -x ) . 

The class of Chebychev-Hermite polynomials are 

2 

orthogonal with respect to exp(-x ) over the real line, and 
are defined by [5] 

^ ,n 

Q (x) = (-1)^ exp (x'^) 2— _[exp (-x^) ] (29) 

" dx^ 

The values of x, can be found from the zero's of 

k 

(29) and those of A^ from (24). However these values have 
been tabulated and are included in Table I for values of n 
from 2 to 6 and in [5] for values of n from 2 to 32 . 

The magnitude of the error bound for approximating 
(28) by a Gauss Hermite Quadrature sum can be shown to be 
[5] 

_ n!(TT)^'^^ maxid^^h(x) 

0 

2 (2n) I dx 

3 . Application to the Example 
The integrand for the example, equation (15) is of 

the Gauss Hermite form since the variables all have Gaussian 



(30) 
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TABLE I 



GAUSS HERMITE QUADRATURE 
COEFFICIENTS AND ZEROS 



Nuinber 

of Evaluation 
Points, n 


Zeros 


X . 
1 


Coefficients A. 

3 


2 


+0.707 


1068 


0.886 


2269 


3 


+1.224 

0.0 


745 


0.295 

1.181 


4090 

635 


4 


+1.650 


680 


0.081 


31284 




+0.524 


6476 


0.804 


9141 


5 


+2.020 

+0.958 


183 

5725 


0.019 
0 . 393 


95324 

6193 


6 


+2.350 


605 


0.004 


530010 




+1.335 


849 


0.157 


0673 




+0.436 


0774 


0.724 


6296 
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TABLE II 



GAUSS LEGENDRE QUADRATURE 
COEFFICIENTS AND ZEROS 



Number 

of Evaluation 
Points, n 


Zeros 


X . 
1 


Coefficients A. 

3 


2 


+0.577 


3503 


1.000 


0000 


3 


+0.774 

“o.o 


5967 


0.555 

0.888 


5556 

8889 


4 


+0 .861 


1363 


0.347 


8548 




+0.339 


9810 


0.652 


1452 


5 


+0.906 


1798 


0.236 


9269 




+0.538 


4693 


0.478 


6287 




o 

• 

o 




0.568 


8889 


6 


+0.932 


4695 


0.171 


3245 




+0.661 


2094 


0.360 


7616 




+0.238 


6192 


0.467 


9139 
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distributions. By changing variables, let 



X 



2 



dx' 



dx^ 





/2o^ 




( 31 ) 



(32) 



(33) 



(34) 



Since X 2 and are independent, the double integral 
becomes the product of the two integrals and (15) takes the 
form 



.00 



#00 

2 2 

g (z , x^ / xp exp ( -X2 -x^ ) dx^dx^ 

— 00 



n n 

= I I A.A, g(z,x ,x ) 
j=l k=l j h 



where 



(35) 



g(z,x^,x^) 



exp [-{z ( 2 ^'^^ 02 X 2 +X 2 ) -1/z {2^^^ o -x^}^ / 2a ^ ] 



(36) 



The value of n can then be chosen as desired and the sum 

and X- and from 

k -^k 

Table II or [5]. The expression (36) becomes quite formidable 



in (35) formed using the values of Aj^ 
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to differentiate after two or more differentiations, so that 
one must assume that the error bound (30) is small. The 
approximation (35) has been performed for n = 6 and the mean 
values R = 1000 ohms, L = 1.0 millihenry, C = 0.001 micro- 
farad. The standard deviation of each variable was assumed 
to be 3.33% of its mean value. The results are contained in 
Chapter IV. 

Not all variables have Gaussian distributions, and most 
performance measures, if functions of many variables, are 
implicit functions. If the example performance measure is 
considered to have been obtained implicitly, with the same 
density functions for the components, the performance meas- 
ure's cumulative distribution function can be determined by 
approximating (19) by a Gauss Legendre Quadrature sum. Since 
exp (-4. 5) is very nearly zero, the limits of and “> can be 
changed to -3a to 3a for each variable and by the following 
change of variables, normalized to the interval [-1, 1]. Let 





i = 1, 2, 3 



(37) 



dx! = dx./3a. 

X 11 



(38) 



so that (18) becomes 



F(z) = I 



Z <z 




(39) 



-1 -1 



where 
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g(xj^,x^ ,x^) 



2 7exp ) 



(40) 




The approximation of (39) by the summation 



Z <^z 




(41) 




evaluating g at these values of x| . It must be remembered 



otherwise g(x') = 0. 

Equation (41) is evaluated for n = 6 and the same values 
of the components as above for the Gauss Hermite example, and 
the results are. contained in Chapter IV for comparison. As in 

j 

the Gauss Hermite example, the error is assumed to be small 
due to the factorials in the denominator of (27) . 



T 

that g takes the form of (40) only for Z (x ) _< z for some z; 
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III. MONTE CARLO METHODS 



This chapter examines the most crude of the Monte Carlo 
techniques for obtaining statistics on a performance measure, 
to methods which attempt to approximate single integrals in 
such a manner that the approximation is reasonably accurate. 
The approximations are then extended to apply to multiple 
integrals in such a manner that fewer evaluations of the 
integrand are made than in conventional methods. 

Monte Carlo methods consist of solving problems of a 
computational nature by constructing a random process for the 
problem, such that the parameters of the random process are 
the desired quantities of the problem. 

Random processes usually imply random variables which 
must be generated in a manner as to represent typical values 
from the variable's probability density function. Random 
number generation is then an important aspect in any Monte 
Carlo method of approximation. Although random number gener- 
ation is a field of its own. Appendix A treats the subject 
suitably for purposes of this thesis. 

A. CRUDE MONTE CARLO 

The most crude form of Monte Carlo is that of taking sam- 
ples of the performance measure to obtain an expected value 
and perhaps a variance of the performance measure. If the 
component-variables are generated in such a way so that the 
values are representative of their distribution functions, and 



24 



for each set of values of component variables generated, 

the performance measure = Z(X^) is evaluated, then the 

expected or average value of the performance measure is [7] 

z = E(z.) = i j 2 (42) 

1=1 

where n is the number of samples. 

While for sufficiently large n, Z is an accurate estimate 
of the average value of the performance measure, no knowledge 
of the form of the probability density function of Z is de- 
rived; and consequently no knowledge of the likelihood of 
other values of Z is obtained. 

A more useful scheme would be to order the sample values 

of Z . as follows to obtain a cumulative distribution for the 
1 

performance measure. 

Consider a random variable S such that if a value of the 
performance measure Z is sampled, and the sample value Z^ is 
less than an arbitrary value, say a, then = 1; and if Z^ 
is greater than the value a, = 0. Then S represents a 

T 

random process. If the values of the component-variables X 
are sampled from their density function as above, then the 
cumulative distribution function for Z can be approximated by 

1 ^ 

F(a) = P[Z^a] = - Is. (43) 

^ i=l ^ 

where n is the number of times the performance measure is 
sampled. 

Since each of the component-variables are sampled inde- 
pendently for each sample value Z^, the number of samples is 
not directly dependent on the number of variables. 
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This crude method yields a cumulative distribution func- 
tion which may be fitted to a polynomial and differentiated 
analytically to obtain the probability density function for 
the performance measure Z. 

As an example of the two methods, consider the example 
of Chapter II. Since the three variables X^, X 2 , and X^ are 
assumed to be Gaussian, three variables w^^, W 2 , and w^ are 
generated by the method in Appendix A to approximate Gaussian 
distributions, each with 0 mean and variance of 1. The val- 
ues of the samples x^, x^ are then 



X. = w.a. + X. , i = 1,2,3 
1 XI X ' ' 



(44) 



where and x^ are the standard deviation and expected val- 
ue, respectively of the variable w.. These values are then 
used to evaluate the performance measure. If the expected or 
average value of the performance measure is desired, the cal- 
culated values of the performance are averaged as in (42). 

If a cumulative distribution is desired, then for k values of 
, j=l,2,...,k the random variable S defined above is summed 

to form F(a.). These values for the examole are listed in 
3 

Chapter IV, and a computer program listing for the computa- 
tion is provided in Appendix B. 



B. CONFIDENCE INTERVALS 

The crude Monte Carlo method above, as with all Monte 
Carlo methods, yields results which are as accurate as de- 
sired within a probability or confidence determined by the 
desired accuracy, and the method of approximation used. What 
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follows is a statistical estimate of the error introduced by 
use of Monte Carlo approximations. 

If A is an event (such as the event that the values of a 
performance measure is less than some specified value) , and R 
is a random variable, such that = 1 if A occurs on the ith 
sample, and R^ = 0 otherwise, and if T is a random variable 
such that 

n 

T = ^ R. (45) 

i=l ^ 

for n samples, then T is the number of times that the event A 
occurs in n samples. The expected value of T is 

n 

E(T) = E( ^ R. ) = nE(Rj = np (46) 

i=l ^ ^ 

where p is the probability of the event A occurring. The 
variance of T is 

n 9 9 

V(T) = V( R. ) = n a (47) 

i=l ^ 

2 

where a is the variance of the event A. 

From Chebychev's inequality [7], 

P = P[lT-npl< d] > 1 - o^/nd^ (48) 

where d is an arbitrarily small number. Equation (48) states 
that for a fixed confidence of [T-np[ being less than a small 
value d, the number of samples to achieve the accuracy, d, 
is bounded by 

2 

n > (49) 

d^(l-P) 

where P is the desired confidence or probability. 
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mum 






Another way of looking at (48) is that for a fixed con- 

■ - 1/2 
fidence P, the error of the estimation varies as on . If 

the variance of the random variable can be reduced to 0, the 

confidence of the estimation could be 1 for an error of 0. 

For crude Monte Carlo, the random variable S, which takes 

on values of 1 whenever the sampled performance measure is 

less then a specified value, and values of 0 otherwise, has a 

binomial distribution. If p is the expected value of S, then 

2 

the variance of S is a = p(l-p) [7]. Thus for a confidence 
of P = 0.9 of being within d=0.01 of np for crude Monte Carlo, 

^ - (0.0001^(0.1) "" p(l-p)10 (50) 

Since the maximum value of p(l-p) is 0.25, inequality (50) 
states that n must be greater than twenty-five thousand in 
order to insure the accuracy of 0.01 with a confidence of 

0.9. 



C. VARIANCE REDUCTION BY ESTIMATING INTEGRALS 

As was seen in Chapter I, the difficulty of determining a 
performance measure's probability density function lies in 
evaluating the integral (1) . Two routines were presented 
which approximated the integral with some accuracy. The 
crude Monte Carlo method illustrated above can be thought of 
as an approximation to an integral. Other methods exist, 
however which reduce the variance by hundreds of thousands 
over crude Monte Carlo [3] . These methods are presented below 
for single integrals and then a general extension is made for 
multiple integrals. 
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If I is the estimator for the integral, f(x)dx, where 

J 0 

X is a variable normalized to the interval [0, 1], then by 
letting 



T _ 

Py(y) 



(51) 



where p^Cy) is the probability density function of the ran- 
dom variable Y and f(y) is the value of the function f at 
some point y sampled from PY(y ) r the expected value of I is 



E(I) 



= E[ 



f (y) 






f (y) dy 



(52) 



It is then necessary to select the probability density func- 
tion Py(y) such a manner as to minimize the variance. 

Define a function f * (x) such that the function which is 
f(x)-f*(x) is a straight line passing through the end- 
points of f(x) on the interval [0, 1], as shown in Figure 2. 
That is, 

f*(x) = f(x) - (l-x)f(O) - xf(l) (53) 

The difference function f(x)-f*(x) can be readily inte- 
grated and can be used as a first approximation to 
rl 

f (x) dx as 

;1 

[f (x) -f*(x) ]dx = l/2f(0) + l/2f(l) (54) 

•• 0 



Now if I* is defined as 



I* = 



(55) 



it becomes necessary to sample f*(y) to get the estimate from 
(53) , (54) and (55) for I 
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f 




Figure 2. Difference Curve Estimation. 



30 



I 



(56) 



= f(0) + f(l) . f (y) - (1-y) f (0) - y(l) 

2 Py(Y) 

To determine p^Cy) / let f(x) be quadratic in x, that is 

f (x) = Ax^ + Bx + C (57) 

By substituting for f in (56) the f of (57) , and replacing I 
with the integral of f in (57) , that is 

I = A/3 + B/2 + C (5 8) 

one can obtain for Py(y) / 

= 6y(l-y) (59) 

It can be sho^m by substitution that the approximation I in 
(56), with Py(y) of (59), is a zero variance estimator for 
the integral of quadratic integrands. 

The function f * (x) could have been defined as 



f*(y) = f(l-y) - yf(0) - (l-y)f(l) 



with a resulting estimator of 

^ f(0) + f(l) . f(l-y) - yf(0) - 

2 2 6y(l-y) 



(60) 



(1-y) f (1) 



(61) 



which is also exact for f (x) quadratic in x. Since y and 
1-y have the same distribution, but lie on the opposite side 
of 1/2, the average of (56) and (61), as can be shown by 
substitution, results in a zero variance estimator for f(x) 
cubic in x[3]. Thus 

^ _ f(0) + f(l) . f(y) + f(l-y) - f(0) - f(l) 

^3 " 2 r2yTP'7) ^ ^ 

is exact for f(x) cubic in x. 
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I 




In the same manner as above, by defining f* as 

f*(x) = f (x) -d-x) (l-2x) f (0)+x(l-2x) f (1) -4x(l-x) f (i) 

(63) 



a difference curve which is quadratic and passes through the 
two end points and the center of f(x) becomes the basis for 
the initial estimate for the integral. Applying the methods 
above, a zero variance estimator for quartic integrands be- 
comes 13]. 

, _ f (0)+f (l)+4f (1/2) 

4 6 



+ f (y) -d-y) (l-2y) f (0)+y(l-2y) f (1) -4y(l-y) f (1/2) 

30y(l-y) (l-2y)^ 



(64) 



and a zero variance quintic estimator becomes 2 

0 

T _ f (0)+f (l)+4f (1/2) 



+ f (y)+f (1-y) -d-2y) ff (0)+f (Dl-Sy (l-y) f (l/2) 

30y(l-y) (l-2y)^ 



(65) 



The probability density function for both the quadratic 
and cubic zero variance estimators is (59) while that for the 
quartic and quintic zero variance estimators is 

Py(y) = 30y (l-y) (l-2y)^ (66) 



D. EXTENSIONS TO MULTIPLE INTEGRALS 

The zero variance estimators above were derived in [2] 
only for single integrals. I'Thile the method can be extended 
to multiple integrals by nesting, nothing is gained since the 
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number of function evaluations increases exponentially with 
the number of variables of integration. In this section, 
double integral extensions to the quadratic, cubic, quartic, 
and quintic zero-variance estimators, and triple integral 
extensions to the quadratic and cubic estimators are proposed, 
and the method of derivation is presented. Although not 
thoroughly tested on a wide range of integrands, it is be- 
lieved that this method of extension can, for multiple inte- 
grals with many variables of integration, provide an accurate 
means of approximating the integral with fewer evaluations of 
the integrand than for the quadrature sum routines. 



For the double integral 




f (x,y) dxdy. 



consider the 



difference function which is a surface passing through the 
four corners f{0, 0), f (0 , 1), f(l, 0), and f(l, 1). This 
function is 



f (x,y) -f*(x,y) = (1-x) (1-y) f(0,0) + (l-x)yf (0,1) 

+ x{l-y)f(l,0) + xyf(l,l) (67) 

Integrating this difference function results in 






[f (x,y) -f* (x,y) ] dxdy = 



f (0,0)+f (1.0)+f (0,l)+f(l,l) 



( 68 ) 



Now if I* is defined similar to (55) as 



I* = 



f*(x,y) 

Pj,(x)p^(y) 



( 69 ) 



and p^(x) and p^Cy) have the form of (59), solving for 
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f*(x,y) in (67) and using (69) as the estimate for 



•1 rl 



•1 fl 



•' 0^0 
fl fl 



•• f * (x,y) dxdy , the estimate to f(x,y)dxdy becomes 



0 J 0 



f(x,y)dxdy = 



f (x,y) - (1-x) (1-y) f (0 ,0) - (1-x) yf (0 , 1) -x(l-y) f (1 ,0) -xyf (1,1) 



36xy(l-x) (1-y) 



(70) 



which is the double integral extension to (56). Replacing x 
with 1-x and y with 1-y in (70) and averaging the resulting 
expression with (70) gives the double integral extension to 
the cubic zero variance estimator (62) 



>1 

0 



f- 



.r, ^ f (0,0)+f (0,l)+f (l,0)+f (1,1) 

f(x,y)dxdy = — ^ — - — ^ — — - 



, f (x,y) +f (1-x, 1-y) - ( l-x-y+2xy) [f(0,0)+f(l,l) ] 

72xy(l-x) (1-y) 

-(x+y-2xy) [ f (0 , 1) +f ( 1, 0 ) ] 

72xy(l-x) (1-y) 



For higher order extensions, a difference surface is 
defined such that it passes through the 2 end points of 
f (x ) , where x is a k vector. This difference curve, which 
is a polynomial in each of its variables, can be integrated 
analytically to give a first approximation to the desired 
integral. The estimate is then refined by sampling the f* 
function weighted by the distribution of the variables, as in 
(69). Thus the triple extension to the quadratic estimator 
(56) becomes 
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f(x,y,z)dxdydz = < (° ■ 0 . 0) (0 . 0 , 1) (0 , 1 , 0) +f (0 , 1 , 1) 

8 



■ 0 



-1 

0 



+ f (1> 0,0) +f( 1,0,1) +f(l, 1,0) +f (1,1,1) ^ 



f (x,y , z) 

216xyz (1-x) (1-y) (1-z) 



, -d-x) (1-y) (1-z) f (0,0,0) - (1-x) (1-y) zf( 0,0,1) 

216xyz(l-x) (1-y) (1-z) 



-(l-x)y (1-z) f (0,1,0) 
216xyz (1-x) (1-y) (1-z) 



, -( 1-x) yz f (0,1,1) -x( 1-y) (1-z) f ( 1 , 0 , 0 ) -x (1-y ) z f ( 1 , 0 , 1 ) 

216xyz(l-x) (1-y) (1-z) 



^ -xy (1-z) f (1,1,0) -xyzf (1,1,1) 

216xyz(l-x) (1-y) (1-z) ^ 

and the triple extension to the cubic estimator (62) becomes 



f(0,0,0)+f(0,0,l)+f(0,l,0)+f(0,l,l) 

8 



4 . f (1, 0,0) +f(l, 0,1) +f(l,l\0 )+£(!, 1,1 ) , f (x,y,z)+f (l-x,l-y,l-z) 

8 432xyz(l-x) (1-y) (1-z) 



- ( 1-x-y-z+xz+yz ) [f(0,0,0)+f(l,l,l)] 
432xyz(l-x) (1-y) (1-z) 



- (z-xz-yz+xy) [ f ( 0 , 0 , 1) +f ( 1 , 1 , 0 ) ] 
432xyz(l-x) (1-y) (1-z) 



- (y-xy-zy+xz) [f(0,l,0)+f(l,0,l) ] 
432xyz(l-x) (1-y) (1-z) 



- (x-xz-xy+yz) [f(l,0,0)+f(0,l,l) ] 
432xyz(l-x) (1-y) (1-z) 



For the quartic and quintic double integral extensions, 
the difference surface must additionally pass through the 
four points which are peripheral midpoints and the internal 
midpoint of the integrand. The double integral extension 
for the quartic zero variance estimator becomes 



r ii 

J n J 0 J 



f (x, y , z) dxdydz 
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rl rl f (0,0)+f (0,l)+f (l,0)+f (1,1) 

f(x,y)dxdy = 

0 •' 0 36 

f (0,y)+f (l,i)+f (^,0)+f (i l)+4f (y,y) 

+ 4 f ± ± ± f ■] 

36 



f (x,y) -(1-x) (l-2x) [ (1-y) (l-2y) f (0,0) 
900xy(l-x) (1-y) (l-2x) ^ (l-2y) ^ 



+ 4y(l-y) f (0,j) -y(l-2y) f (0,1) ] 

900xy(l-x) (1-y) ( l-2x) ^ ( l-2y) ^ 

-4x(l-x) [ (1-y) (l-2y) f (|,0)+4y(l-y) f (i,|) -y(l-2y) fj,l) ] 
900xy(l-x) (1-y) ( l-2x) ^ ( l-2y) ^ 



x(l-2x) [ (1-y) (l-2y) f ( 1 , 0 ) +4y ( 1-y) f ( 1 , y) -y ( l-2y ) f ( 1 , 1) ] 

+ ^ ^ 

900xy(l-x) (1-y) (l-2x) ^ ( l-2y) ^ 

and the double integral extension for the quintic zero 

variance estimator becomes 






f (0,0)+f (0,l)+f (l,0)+f (1,1) 
36 



f (0,y)+f (l,y)+f(y,0)+f (y,l)+4f (y,y) 

+ 4[ f £ f £ ±^J—] 

36 



f (x,y)+f (l-x,l-y) -d-2x) (l-2y) (l-x-y+2xy) [ f (0 , 0 ) +f ( 1 , 1) ] 

1800xy(l-x) (1-y) ( l-2x) ^ (l-2y) ^ 

(l-2x) (l-2y) (x+y-2xy) [ f (1,0) +f (0 , 1) ] 

1800xy(l-x) (1-y) ( l-2x) ^ ( l-2y) ^ 

-4y(l-2x) (1-y) [f (0,j) ^ 

1800xy(l-x) (1-y) (l-2x) ^ (l-2y) ^ 

-4x(l-x) (l-2y) [f (y,0) -f (y,l) ] -32xy(l-x) (1-y) f (y,y) 

+ ^ 5 5 (75) 

1800xy(l-x) (1-y) (l-2x) ^ (l-2y) ^ 
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For higher order extensions to the quartic and quintic 

estimators, the difference function must be defined so as to 

}c }c “ X 

pass through the 2 corners, the k.2 peripheral midpoints, 

T T 

and the central midpoint of f(x ) , where x is a k vector. 
These extensions become difficult to derive due to the large 
number of points involved. 

The triple integral extension to the cubic estimator, (73) 
has been applied to the example of Chapter II. The results 
are contained in Chapter IV for comparison, and a printout 
of the program is contained in Appendix B. 
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IV. CONCLUSIONS 



Chapter II presents the mathematics for obtaining a 
performance measure's probability density function for ex- 
plicit and implicit performance measures. Two integration 
schemes using Gauss quadrature sums were then discussed. 
Chapter III illustrated a crude Monte Carlo method for ob- 
taining a performance measure's distribution cind then dis- 
cussed the confidence of such schemes. Monte Carlo schemes 
for evaluating single integrals were discussed. Double and 
triple integral extensions to these methods were then derived. 

From observation of the problem, it seems intuitive that 
most system's performance measures will be of the implicit 
type, and in general the performance measure will be a func- 
tion of several variables. In the process of evaluating the 
performance measure's distribution, the integral (1) must be 
evaluated. Evaluation by the quadrature sum methods, (25) 
and (29) require some n evaluations of the integrand, where 
n is the number of points in the quadrature sum, and k is the 
number of variables. This number can become quite large, and 
since the computation time for such routines is roughly pro- 
portional to the number of times the integrand is evaluated, 

« 

this process could take a good deal of time. However, this 
method is quite accurate, and it is easily programmed for a 
large number of variables. 
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The reduced variance estimators reduce the number of 



evaluations of the integrand required. For a single esti- 

mate, the quintic extension requires evaluations at the 2 

k-1 

corners, the k2 peripheral midpoints and the 3 internal 

k-1 

points, or a total of 3 +(k+2)2 points, where again k is 

the number of variables. For additional estimates, only two 

additional points per estimate are required. Thus M estimates 

k -1 

of the integrand would require 2M + 1 + (k+2)2 evaluations 
of the integrand. This can, in most cases be much less than 
that required for the quadrature sums. 

The drawback to the reduced variance scheme is that the 
geometric form of the integrand must be examined to insure 
that the major portion of the surface will not be neglected in 
the sampling process. As illustrated in Figure 3, the dif- 
ference curve does not provide a good representation of the 
function, while splitting the curve into two sections as in 
Figure 4 enables one to obtain a good representation of the 
function by the use of tv/o difference curves. 

The generation of the random numbers with the probability 

2 

density function 30x(l-x) (l-2x) can be bothersome to program 
as can the quintic estimator extension. The cubic estimator 
is much more easy to extend, and the density function 6x(l-x) 
is more readily programmable. However more points are re- 
quired by the cubic to ensure as good accuracy as the quintic 
extension . 

Properly used, the Monte Carlo method of integral estima- 
tion can provide a more rapid method of determining a 
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Figure 3. Improper Difference Curve Choice. 
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Figure 4. Proper Difference Curve Choice. 
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performance measure's density function with reasonable ac- 
curacy than other methods which require more standard inte- 
gration routines. 

A. RECOMMENDATIONS FOR FURTHER STUDY 

Monte Carlo schemes seem to be attractive for reducing 
the amount of time required in the evaluation of the multi- 
ple integrals required in evaluation of a performance 
measure's distribution. Higher order zero variance estima- 
tors could be generated by fitting a difference curve to the 
points which are the zero's of a Legendre polynomial and 
evaluating the required probability density function to 
achieve a zero variance estimator for integrands of order 
2n-l, where n is the order of the Legendre polynomial. Other 
schemes of generating random numbers of particular distribu- 
tions could as well be investigated. 

B. RESULTS OF APPLICATIONS TO EXAMPLE 

The cumulative distribution function of the upper 3 db 
frequency of the series R-L-C circuit of Figure 1 has been 
obtained by four methods. 

First, using the analytical expressions (11) and (12) for 
the integral (6), the probability density function, equation 
(15) was obtained by evaluating the integral with a 6 -point 
Gauss Hermite quadrature sum routine. The cumulative dis- 
tribution function was then obtained from the probability 
density fixnction by rectangular integration. The resulting 
distribution has been plotted in Figure 5, and will serve as 
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the reference since integration of explicit functions is 
generally more accurate than integration of implicit 
functions . 

Using the same performance measure as for the explicit 
case, but stipulating that the value of the performance meas- 
ure was to be obtained implicitly, integral (18) was evalu- 
ated using 6 -point Gauss Legendre quadrature sums. As can be 
seen in the plot of the cumulative distribution function for 
this case. Figure 6, the distribution is very nearly the same 
as that obtained in the explicit case. 

Then for the crude Monte Carlo case, 1000 samples of the 
performance measure (10) were taken with the values of the 
three variables being obtained from a Gaussian random number 
generator. As can be seen in Figure 7, the resulting cumula- 
tive distribution function is very nearly the same as that 
obtained in the explicit case. 

The extended cubic zero variance estimator equation (73) 
was then applied to the integral (18) . The required distri- 
butions of the components were generated as explained in 
Appendix A. Twenty approximations for each value of the per- 
formance measure averaged to obtain the cumulative distribu- 
tion plotted in Figure 8. Two discrepancies are observed 
here. First, the cumulative distribution exceeds the maximum 
value of 1, and secondly, there is an apparent discontinuity 
in the curve. Since the entire curve has a higher value than 
that of Figure 5, it is assumed that the difference surface 
estimation was too large. 
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Angular Frequency xlO® 

Figure 5. Cumulative Distribution Function for Example Using Gauss Hermite Quadrature 
(Explicit) . 



Probability 
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Angular Frequency xlO^ 

Figure 6. Cumulative Distribution Function for Example Using Gauss Legendre Quadrature 
(Implicit) . 
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Angular Frequency xlO® 

Figure 7. Cumulative Distribution Function for Example Using Crude Monte Carlo. 
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Although the example did not prove the worth of the 
Monte Carlo method, a count of the required number of evalu- 
ations for high multiple integrands for a 6 -point Gauss 
Legendre routine as compared to that required in the extended 
cubic estimator averaging 100 approximations for each value 
of the performance measure will illustrate its merits. For 
a multiple integral with 10 variables of integration, the 
6-point Gauss Legendre method would require 6^^ evaluations 
of the integrand. For the same integrand, the cubic exten- 
sion would require 20 0 + 2^*^ evaluations. The savings in 
computation time should be apparent. 
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Appendix A 



RANDOM NUMBER GENERATION 
A. UNIFORM RANDOM NUMBERS 

A truly uniform random number cannot be generated on a 
digital computer. However, methods such as the power residue 
[4, 10] method can generate a sequence of numbers which have 
the properties of uniformly distributed random numbers. 

The method of power residues is based on modular arith- 
metic. If is an arbitrary integer in the sequence of 
pseudo-random numbers to be generated, the next number gener- 
ated is determined by 

= X^p (modulo M) A-1) 

2 

where p is any prime integer such that p is greater than M, 

the modular base of the computer's arithmetic unit. The 

mamber X^/M is then an approximation to a random variable 

with a uniform distribution on the interval (0,1). Utilizing 

a digital computer with a 32 bit register (31 number bits 

plus one sign bit) , a value of p of 65,539 for an M of 

31 29 

2,147,483,647 = 2 -1, a sequence of 2 numbers can be gener- 

ated before any number of the sequence is repeated [4] . 

Since the digital computer automatically performs modular 
arithmetic on integers, with the modulus being determined by 
the size of the registers in the arithmetic unit, uniform 
pseudo random numbers can be generated quite rapidly. 
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B. GAUSSIAN RANDOM NUMBERS 



Gaussian random numbers are most commonly generated uti- 
lizing uniformly distributed pseudo-random numbers. If 
R^ , i=l, 2 ,...,n are independent samples from a density func- 
tion uniformly distributed over the interval ( 0 , 1 ), then a 
nomally distributed random variable with 0 mean and vari- 
ance of 12 /n can be approximated by [3, 4, 10] 

1 0 ^ 

2 = ~ I (R.) - 6 (A-2) 

i=l ^ 

If n = 12, Z will have 0 mean and variance of 1. 

C. 6 x(l-x) DISTRIBUTION 

If R^ is an independent sample from a random variable 
uniformly distributed on the interval (0,1), for i = 1,2,3, 
and if the three samples are rearranged in order of magni- 
tude such that 

Rj^ <. R 2 - ^3 (A-3) 

then R^ has probability density function 6 x(l-x) [3]. 

D. 30x(l-x) (l-2x)^ DISTRIBUTION 

If R^, R 2 , 1 ^ 2 ' ^ 4 ' ^5 independent samples from a 

random variable uniformly distributed over the interval ( 0 , 1 ) , 
and are arranged so that 



^1 2 1 - 1^2 2 1 - * 



< IR5 - j\ 



(A-4) 



and if S is chosen such that S = R^ with probability 3/4 and 
S = R^ with probability 1/4, then S has probability density 



50 



I 



2 

function 30x(l-x) (l-2x) [3], This density function, it was 

discovered, requires some complex programming to generate. 



E. OTHER DISTRIBUTIONS 

Let Y be a uniformly distributed random variable over the 
interval (0,1), and let 



where h(t) 



y = F(x) 

1 for 0 t 1, 

X 

F(x) = f 



= j h(t)dt 

..oo 

and h(t) = 0 

Pj^(x) 



otherwise . 



(A-5) 

Since 

(A-6) 



where p^Cx) is the probability density function for X, if 
F(x) or p (x) , the desired function, are known explicitly, 
then the inverse function 



X = F"^(y) 



(A-7) 



will generate x according to its desired density function 
p^(x) [3]. Most often p^ (x) will be a polynomial fitted to 

some statistical data which represents the distribution of 
the variable X. Since this polynomial is normally of a high 
order, some difficulty can result from this method in that 
the zeros of (A-6) may be difficult to find. 
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<^PPEMOIX B 



COMPUTER PROGRAMS 



COMPUTATTOM CP CUMULATIVE DISTPIBUTICN pUNCTIOM OF 3 
CB AMGULAP F^EQUEMCY OF THE SERIES R-L-C CIRCUIT 
EXAMPLE. P'=R=rcMAMCE 'lEASUoP IS AN EXPLICIT puNCTION 
OP THE VART^BLES. ‘^-THOO IS 6 POINT GAUSS HPRMITE 
CUACRATURE SUM APPROXIMATION TO EQUATION (15» TO 
OBTAIN FPCBABILITY CENSITY FUNCTION. RECTANGULAR 
INTEGRATION IS THP'N USED TO OBTAIN CUMULATIVE 
DISTRIBUTIOk F’imcTICN. 

COMMON Z 

DIMENSION a(6) ,Y(6) ,X( 2) 

WRITE(6 .61 ) 

F = 0. 

INPUT GAUSS HER'MTE COEFFICIENTS AND ZEROS 
A(l)=.A53001E-2 
A( 2)=. 15706C7 
A(3)=.72A6296 
A (4) =A(3) 

A ( 5 ) =A ( 2 ) 

A(6)=A(1) 

Y(l) =2.3506C5 
Y( 2)= 1. 335849 
Y (3)=.43607A1 
Y(4)=-Y(3) 

Y(5)=-Y( 2) 

Y(6 l=-Y ( 1 ) 

INPUT WCPST CASE FREOUFMCIcS AND COMPUTE DELTA Z AND 
CISCRETE FCPOUENCIES. 

ZHI=. 18791S5E7 
ZLO=. 14C5987E7 
PELZ=(ZHI-ZL0)/19. 

Z = ZLO 

FOR '^ACH Z, EVALUATE PDF USING 6 POINT GAUSS HERVITE 
CUADRATUPP . 

DC 3 I Z=1 ,2C 
P=0. 

DC 2 1=1,6 
DO 2 J=1 ,6 
X( 1 ) = Y( J ) 

X(2 )=Y( I ) 

2 P=P+A( I )-':‘A( J)*=UN( X) 

USE rectangular integration now to obtain cumulative 

DISTRIPIJTICN. 

F = F + P'-.'!D£LZ 
C OUTPUT CU''ULATIVF DISTRIBUTION 

WRITE (6 .60) ZjF 
C INCREMENT Z 

3 Z=Z+DELZ 
STOP 

60 PORMATI 2F30. 7) 

61 FCRMATI ’ 1’ .////////, 15X. ’ANGULAR FREOUS NC Y • . 1 6X . 
I'CISTRISUTICN* ,//) 

END 
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PUNCTinN FUN'(X) 

C FlJNCTinM SUBPQIJTIMk for US£ with gauss HERVITE 

C QUADRATURE SUM ROUTINF. 

CGMMON Z 

DIMENSION X ( 2 ) , XO ( 3 ) , S GMA ( 3 ) , X 1 ( 3 ) 

CUN=0. 

C CHECK flag to see T c CONSTANTS HAVE BEEN COMPUTED. 

IF( 11-121 )1 ,3, 1 

1 PI=SORT( 3. 14159) 

P2=SORT(2. ) 

C INPUT MEAN VALUES C= COMPONENTS 

XC( l) = l.R-3 
XO( 2 )=1 .E-9 
X0(3)=l.?3 

C CALCULATE STANDARD DEVIATIONS OF COMPONENTS 

DO 2 1=1.3 

2 SGMAd ) = XC(T)/30. 

C SET PLAG 

I 1 = 121 

3 CGNTINUE 
DO 4 1 = 1 ,2 

C "UMNORMal IZE" VARIABLES FOP. COMPUTATION OF DERIVATIVE 

C AND PiJNCTICN INVEPSE. 

4 X 1 ( I ) = X( I ) *p 2=;= SG '•’A ( I » + XO ( I ) 

X1(3J = Z+-X1( 1 )-!./( Z*X1( 2) ) 

DRV=X1 ( 1 ) + l . / ( Z=i= Z^=X1 (2 ) ) 

C COMPUTE FX?nM5MT 0= ^'TI AL 

F=(X1(3)-X0(3) )/ (P2^SGMA( 3 ) ) 

F = c * c 

IP( F-13.6) 5,6, 6 
C EVALUAT- FUNCTTOM 

5 FIJN = DRV^'=L-XP(-F)/(P2'^^I*SGMA(3) ) 

6 RETURN 
END 
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GULAR FREOL 


IcMCY 


01 STRIBUTION 


0. 1405C87.1 


C7 


0.93585405-06 


0 . 1A3C892E 


07 


0.25853635-04 


0.14557S7E 


07 


0 .20947625-03 


0. 148C7C2F 


C7 


0. 23163275-02 


0.15056C7E 


07 


0. 13032755-01 


0.1530512= 


07 


0.49607905-01 


0.1555417E 


07 


0.1364721E 00 


0.15803225 


07 


0.29398995 00 


0.1605227= 


07 


0.49353155 00 


0. 163C132E 


07 


0.6891426= 00 


0.1655037E 


07 


0.84302755 00 


0. 16799425 


07 


0.9314889= 00 


0.17048475 


07 


0.97439035 00 


0.17297525 


07 


0.99251915 00 


0.17546575 


07 


0.99806685 00 


0.1779562E 


07 


0.9994828F 00 


0. 18C44675 


07 


0.99988135 00 


0.18293725 


07 


0.9999751= CO 


0.18542775 


07 


0.9999889= CO 


0. 1879182E 


07 


0.9999900= 00 
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CCVPUTATICN OF CU'^ULATIVE DISTRIRUTION FUN'CTIGN OF 3 
CB ANGULAR =REOU£NCY OF THE SERIES R-L-CLC CIRCUIT 

example, performance measure is implicit function of 

VARIABLES. METHOD IS 6 '^OINT GAUSS LEGENDRE 
QUADRATURE SUM APPROXIMATION TO EQUATION (18). 

COMMON Z(2C),IZ 

CIMENSICN A (6 ) , Y (6 ) ,X(3 ) » F( 20 ) 

WRIT£(6,61 ) 

C INPUT GAUSS LEGENDRE COEFFICIENTS AND ZEROS. 

A(1 )=.17132A5 
A( 2 )=.36C7616 
A (3 )= .4679139 
A (4) =A( 3) 

A( 5 )=A ( 2 ) 

A (6)=A( 1 ) 

Y(l)=. 9324695 
Y(2)=.6612C94 
Y(3)=. 2386192 
Y( 4)=-Y( 3) 

Y(5)=-Y(2) 

Y(6)=-Y( 1) 

C INPUT WORST CASE FREQUENCIES AND COMPUTE DELTA Z AND 

C DISCRETE 'FREQUENCIES. 

ZHI =. 1879185=7 
ZLO=. 14C5987E7 
DELZ=( ZHI-ZLO/19. 

Z(1)=ZL0 
DO 5 IZ=1,20 
IF( IZ.EC.l ) GO TO 1 
Z( IZ)=Z( IZ-U+CELZ 

1 F(IZ)=0. 

C FCR EACH frequency, CCVPUTE 6 POINT GAUSS LEGENDRE SUM 

DO 2 1=1,6 
DC 2 J=l,6 
DC 2 K=l,6 
X( 1 ) = Y( K ) 

X(2)=Y(J) 

X(3)=Y(I) 

A1=A(I )*A(K ) 

2 F(IZ) = F( IZ)+A1=!--FUN(X) 

C OUTPUT DISTRIBUTION 

5 WRITS(6,60)2(IZ),F( IZ ) 

STOP 

6C FORmAT( 2E30.7) 

61 FORMAT! ' 1* .////////, 15X, 'ANGULAR FREQUEMC Y ' . 1 6X , 

1 'DISTRIBUTION' ,//) 

END 
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FUNCTICN FUN’(X) 

C FUNCTICN SUBROUTTN; FOR USE WITH GAUSS LEGENDRE 

C INTEGRATION SCHEME. 

common Z (20) ,I Z 
DIMENSION XC(3) , SGMA( 3) tY( 3) 

DIMENSION X(3) 

PI = SORT( 2.=i=3.14159)**3 
C INPUT MEAN VALUES OF COMPONENTS 

XC(1)=1 .E3 
XC( 2)=1 .r-3 
X0(3) = l .5-9 

C CALCULATE CCMP0N5NT DEVIATIONS 

DO A 1=1,3 
SG‘-'A( I ) = X0( I )/30 . 

C "UNNORMAlI Z'^" VARIABLES 

A Y( I )=3.-SG^'A(I )--)'X( I ) + X0( I ) 

P=Y(1 )/ (2 .*Y(2 > ) 

C COMPUTE VALUE OF FER-CRMANCE MEASURE FOR THESE 

C COMPONENT VALUES. 

ZZ = P + SORT( F=!=P + 1 ./(Y(2)^Y(3) ) ) 

C CHECK TO Sl-F I P5PFCRMANCE MEASURE IS LESS THAN 

C RRFEPENC5 VALUE. 

IF(ZZ.LE.Z(IZ))GC TO 1 
C IF NOT, THE INTEGRAND IS ZERC 

FUN=0. 

GC TO 3 

C 1'= IT IS COMPUTE VALUE OF INTEGRAND 

1 F=0. 

DC 2 1=1 ,3 

2 F=F+X( I ) **2 
FUN=27.*EXP(-4.5*F)/PI 

3 RETURN 
END 
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f-GULAR FREOL 


;ENCY 


DISTRIBUTION 


0. 14C5F87r: 


C7 


0.0 


0 . 1430892E 


07 


0. 10826325-05 


0. 14557S7E 


07 


0.40259775-04 


0. 14807C2E 


C7 


0.61146645-03 


0 . 15056072 


C7 


0.51073995-02 


0.153C512E 


07 


0.2215C09:->01 


0.1555417E 


07 


0. 8334792~-01 


0.15803225 


07 


0.196093C5 00 


0. 1605227E 


07 


0.39814955 00 


0 . 1630 132E 


C7 


0.58928135 00 


0.16550375 


07 


0.79122165 00 


0. 1679S42E 


C7 


0.90197745 00 


0.17048475 


07 


0.93836035 OC 


0. 17297525 


Cl 


0.96852295 00 


0.17546575 


n 


0.98413295 00 


0.17795625 


07 


0.98672635 00 


0. 18044575 


07 


0.9874057= 00 


0.18292725 


C7 


0.9874^385 CO 


0.18542775 


07 


0.98744585 CO 


0. 18791825 


C7 


0.98744595 00 
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CO‘-»PUTATinN OF CUMULATIVE 0 1 SIR I B UTI CN PUNCTICN GF 3 
CB AMGULAR prEOUENCY Qp THE SETIES R-L-C CIRCUIT 

example. method is crude monte carlo, where the 

PERFORMANCE MEASURE IS SAMCLHO lOOC TIMPS. COMPONENT 
VALUES ARE OBTAINED FROM GAUSSIAN RANDOM NUMBER 
GENERATORS. 

CIMENSICN IX(3) 

D1 me NS I ON Z(20) , X (3 ) ,X0(3 ) ,SGMA (3 ) , F( 20 ) 

WRIT£(6, 61 ) 

INPUT V>fORST CASE FREQUENCIES AND COMPUTE DELTA Z 
ZLO=. 14C5S3E7 
ZHI=. 1879185P7 
DELZ=( ZHI-2LC)/ 19. 

Z( 1)=ZL0 

INITIALIZE DISTRIBUTION TO ZERO. 

F(1)=0. 

DO 1 1=2,20 
Z( I )=Z( I-l )+OELZ 

1 F(I )=0. 

INPUT MEAN VALUES OF VARIABLES AND COMPUTE DEVIATION. 
X0(1)=1 .E3 
X0(2 )=1 .E-3 
XC( 3)=1 .E-9 
DO 2 1=1,3 

2 SGMA( I )=X0( I )/30 . 

INPUT KEPMALS FOR RANDOM GENERATICN 
IX(1)=13A95 
IX(2) =54621 
IX( 3)=17A327 

sample PERFCRMAMCE •MEASURE 1000 TI»^ES 
DC 4 1=1 ,1CC0 

GENERATE COMPONENT VALUES FROM GAUSSIAN DISTRIBUTIONS. 
DC 3 J=l,3 

3 CALL GRN( IX( J) , SGM,A(J) ,X0( J) ,X ( J) ) 

CCMPUTE PERFORMANCE MEASURE. 

PP = X(1 )/(2.-X(2 ) ) 

ZZ=PP+SQRT( PP*M2+ li / ( X( 2)*X( 3) ) ) 

FIND DISCRETE VALUES 0^= Z FC'P WHICH PERFORMANCE 
MEASURE IS LESS THAN OR EQUAL TO, AND INCREMENT THE 
DISTRIBUTION COUNTERS. 

DC 4 J=l,20 

IF(ZZ.LL-.Z(J) )F-( J)=P( J) + l. 

4 CONTINUE 

C DIVIDE DISTRIBUTION COUNTER BY 1000 TO OBTAIN 

C DISTRIBUTION 

DO 6 J=l,20 
F{ J)=F( J)=«=l .P-3 
C OLTPUT DISTRIBUTION. 

6 WRITEI6, 6C )Z( J ) , =( J ) 

STCP 

60 FOR»-<AT( 2E3C.7) 

61 FORMAT! • 1« ,////////, 15X, 'ANGULAR FREOUEMC Y* , 16X, 

1 'DISTRI BUT I CM' ,//) 

END 
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SUBRCUTIMt GPN( I Xt S,AM, V) 

C SUeRPUTINS TO GENERATE GAUSSIAN DISTRIBUTED RANDOM 

C VARIABLES. 

A=0. 

DC 1 1=1.12 
CALL UPM(IX,IY,Y) 

IX=IY 
1 A=A+Y 

V=( A-6. )#S+AM 

RETURN 

END 



SUBROUTINE URN ( I X, I Y, YFL ) 

C SUBROUTINE TO GENERATE UNIFORM DISTRIBUTED RANDOM 

C VARIABLES. 

IY=IX=!'6553G 
IF( lY) 1,2,2 

1 I Y=I Y+21474836A7+1 

2 YFL=IY 

YFL=YFL*.4656613E-9 

RETURN 

END 
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I 



ANGULAR FREQUENCY 



0.14C5G80E 07 
0.1430885E 07 
0.14557SCE C7 
0.148C695E 07 
0.15C560CE 07 
0.1530505?: 07 
0.15554103 07 
0.1580315c 07 
0.1605220= 07 
0.16301255 07 
0.1655C30E 07 
0.16799355 07 
0.17C4840E 07 
0.1729-J459 07 
0.1754650E 07 
0.17795559 07 
0.18C4460F: 07 
0.1829365c 07 
0.1854270c 07 
0.18791755 07 



DISTP IBUTION 



0.0 
0.0 
0 . 0 
0.0 

0.69999999-02 
C. 2600000^-01 
0.7799995E-01 
0.21200009 00 
0.38499995 CO 
0.5929999E 00 
0.75499999 00 
0.8839999‘i CO 
0.94899995 00 
0.97899999 00 
0.9949999t5 00 
0.99999999 00 
0.99999995 CO 
0.99999999 00 
0.9999999c 00 
0.9999999E 00 
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C CCMPUTATION np CUMULATIVE DISTRIBUTION FUNCTION OF 3 

C CE 4NGULAR '^REOUENCY 0*= SERIES R-L-C CIRCUIT EXAMPLE. 

C PERFORMANCE MEASU'^S IS AM IMPLICIT FUNCTION OF ITS 

C VARIABLES. METHOD IS TRIPLE INTEGRAL EXTENSICN TO 

C CUBIC ZERO VARIANCE ESTIMATOR APPLIED TO EQUATION (18) 

COMMON Z(20),IZ 

DIMENSION X( 3) . XI (3 ) . Y( 3) . IX( 3) . IXK 3) .F( 20 
V>RITE(6,61) 

C INPUT WORST CASE FREQUENCIES AND COMoyrE DELTA Z AND 

C DISCRETE VALUES OF Z. 

ZHI = . 187E185'"7 
ZLO= .1405987E7 
OELZ=( ZHI-ZL0)/19. 

Z( 1)=ZL0 
CO 1 1=2.20 

1 Z(I )=Zn-l)+05LZ 

C DO FOLLOWING FOR EACH REFERENCE VALUE OF Z 

DO 999 IZ=1.20 

C INPUT i^ERNALS FOR RANDOM NUMBER GENERATION. 

IX( 1 ) = 49A31 
IX(2)=39^163 
IX( 3) =1755381 
SUM = 0 . 

C COMPUTE "CCRNER VALUES" OF INTEGRAND AND GROUP INTO 

C POUR VALUES. 

DC 2 1=1.3 
X(!)=0. 

2 Xl( I )= 1 .-X( I ) 

F1 = FUN( X ) + FUN(Xl ) 

X(3)=l. 

Xl(3)=0. 

F2 = FUN( X ) + FUM(Xl ) 

X( 1)=1. 

Xl( 1 ) = 0 . 

F3=FUN( X) + = LN(X1 ) 

X( 3)=0. 

Xl(3)=l. 

F4=FUN( X) + FUN( XI ) 

C COMPUTE INITIAL ESTIMATE 

AVG=( =l + F2 + =3 + -=4 )/8 . 

C NOV-I GET THE AVERAGE OF 20 REFINEMENTS 

DO 7 JX= 1, 2C 

C GENERATE 3 RANDOM VARIABLES WITH 6X(1-X) DISTRIBUTION. 

DO 800 J=l,3 
DO 6 1=1.3 

CALL URN( IX ( J) , IXI ( J) ,Y ( I ) ) 

6 !X(J)=’X1(J) 

IF(Y(1)-Y(2))100, 1000. 200 

ICO I F(Y( 1)-Y( 3) )700, 1000,500 
200 IF(Y(2)-Y(3) )3CC, 10CC,400 
3C0 IF( Y( 1 )-Y( 3 ) )5?C , 1000, 600 
4C0 X( J ) = Y( 2) 

GO TO 800 
500 X(J)=Y(1) 

GO TO 8C0 
600 X(J)=Y(3) 

GC TO SCO 

7C0 IF(Y(2)-Y(3))4CC,1CCC,600 
800 XI ( J)=l .-X( J ) 

C COMPUTE REFINEMENT 

SUM1=0 . 

SUM1 = SUM1 + FUN(X )+FUN(Xl ) 

SLM1 = SUM1-=1=^( l.-X(l)-X(2)-X(3)+X(l)=^'X(2)+X(l)*X(3) + 
1X(2)+X(3) ) 

SUM1 = SUmi-F2* (X(3 )-X (1 )=!=X(3 )-X (2 )-X (3 )+X ( 1 )=^X(2) ) 
SUM1=S'JM1-P3=.M X(2)-X(1)=!=X(2)-X(3)*X(2) + X(1)=^X(3)) 

SUM1 = SUM1-P4*(X ( 1)-X( 1)*X(2)-X( 1)*X( 3) + X( 2)=:=X(3) ) 

7 SLM = SUM + SUMl/( 432 .=^X{ 1 )=!=X(2)=:--X (3 )*X1 (1) =)^X1 (2 )*X1 (3 ) ) 

C COMPUTE THE VALUE OP CUMULATIVE DISTRIBUTION FOR THIS 

C VALUE CP Z 

P ( IZ ) = . CS’i'SUM+A VG 
F( IZ )=F( IZ )/l .125 
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C OUTPUT DISTRIBUTION 

9S9 VjRITB(6,60)Z( IZ) ,F (IZ) 

1000 STO>^ 

60 FCRVAT(2B30.7) 

61 FORMAT! • 1* ,//////// ,15X, • ANGULAR FR EQUENC Y* , 16X , 
I'CISTRIBUTIOM',//) 

END 



FUNCTION Fijn(X) 

C FUNCTION SUBROUTINE FCR MONTE CARLO METHOD 

COMMON Z( 20 ,I Z 

DIMENSION X0( 3 ) » SG^'' A ( 3 ) » Y ( 3 ) , X( 3 ) 

C CHECK '"LAG TO DETERMINE I'" CONSTANTS HAVE BEEN 

C EVALUATED. 

IF(II-11A)1,3. 1 

1 PI = SQRT ( 2.-3. 14159 )=!'*3 

C INPUT MEAN VALUES OF COMPONENTS 

X0(1 ) = 1 .E3 
XC( 2) = 1 .'-■-3 
X0(3)=l .c-9 

C CQMPUTf" DEVIATIONS OF COMPONENTS. 

DC 2 1=1 ,3 

2 SGMA( n =XC( I ) /30. 

C SET FLAG 

11=114 

3 G = 0. 

C BREAK INTSGCAND INTO 8 FARTS AND SUM! INDIVIDUAL 

C EVALUATIONS. 

DO 6 I = lt2 

C ••UNNCPMALIZ'’-” VAQ TABLES FOR EACH FART. 

Y( !) = (-! )='.‘=’M^3 .^SG’-'AT 1 )*X ( 1 )+X0 ( 1 ) 

DO 6 J= 1, 2 

Y (2 ) = ( -1 )=«=* J«3 .=^SGMA(2)=S^X(2)+X0(2) 

DO 6 K=l,2 

Y(3)=(-l )=i^*'<*3.*SG‘^A( 3) -X( 3) + X0( 3) 

P=Y( 1) / (2.-Y(2 ) ) 

C COMPUTr VALUE 0= PERFORMANCE MEASURE 

ZZ = P + SQPT( + i ./(Y(2)*Y(3) ) ) 

C CHECK VALUE OF ?ER'"0'’M AMCE MEASURE AGAINST REFERENCE. 

1F(ZZ-Z( IZ) )4,4,6 

C IF LESS THAN OR EQUAL TC PEFERENCEt EVALUATE FUNCTION 

4 F = 0. 

DO 5 L=l,3 

5 F==+X(L)**2 
G=G+EXP(-4.5=f=F) 

C IF GREATER THAN REFERENCE. SET FUNCTION ECUAL TO ZERO. 

6 CONTINUE 
FUN = 27. =!=F/FI 
RETURN 

END 



SUBROUTINE URN ( I X. I Y, YFL ) 

C SUBROUTINE FCR GENERATING UNI=DRM RANDOM NUMBERS. 

IY=IX^65539 
IF! IY)1.2.2 

1 IY=I Y+2147483647+1 

2 YFL=IY 

YFL=YFL='=.4656613E-9 

RETURN 

END 
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A? =^P?QU':MCY 


OISTC ^0^' 


0.1405987" 07 


0.89843640 


-07 


0.1^30892' 07 


0.8984364- 


-07 


0.1455797’* 07 


0. 3900373" 


-Oi. 


0.1480702'- 07 


0.4555574" 


-03 


0. 15056C7? 97 


0.2644415" 


-"■2 


0.1550512- O’' 


0.2990560’ 


-Cl 


0.1555917”- 07 


0.1155540“ 


CO 


0. 1580322'' 07 


0. 1728926" 


GO 


0.1605227" 07 


0.2506310- 


00 


0.1630132 1 07 


0.3067204“ 


CO 


0.1655037- 07 


0. 8 ■'46661 


GO 


0.) 6 79042': 


0.9280446" 


CO 


0.1704847" O’? 


0.1000616 ” 


01 


0.1720752" C7 


0. 1047179" 


01 


0.1754 65 7’: 07 


0. 1061230’’ 


01 


0.1779562- 07 


0.1062497- 


01 


0.180466’'l 07 


G. 1062828- 


''1 


0.1829372’’ 07 


0 . 1C62843' 


01 


0.1856277' 0' 


0. 1C62867" 


Cl 


0.18701823 07 


0. 1062367" 


01 
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